home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / SciAn / src / lorenz.f < prev    next >
Text File  |  1994-08-01  |  5KB  |  216 lines
  1. C -------------------------------------------------------------------------
  2. C SUBROUTINE lorenz(x,y,dydx) RETURNS DERIVATIVES dydx AT x FOR THE
  3. C LORENZ ATTRACTOR
  4. C -------------------------------------------------------------------------
  5.       
  6.       subroutine lorenz (x,y,dydx)
  7. C     =================
  8.  
  9.       implicit real*4 (a-h,o-z)
  10.       parameter (nmax=10) 
  11.       dimension y(nmax),dydx(nmax)
  12.  
  13.       common /parm/ sig,b,r
  14.  
  15.       dydx(1)=sig*(-y(1)+y(2))
  16.       dydx(2)=-y(1)*y(3)+r*y(1)-y(2)
  17.       dydx(3)=y(1)*y(2)-b*y(3)
  18.  
  19.       return
  20.       end
  21.  
  22.  
  23. C -------------------------------------------------------------------------
  24. C GIVEN VALUES FOR n VARIABLES y AND THEIR DERIVATIVES dydx KNOWN AT x, USE 
  25. C THE FOURTH ORDER RUNGE-KUTTA METHOD TO ADVANCE TO THE SOLUTION OVER AN
  26. C INTERVAL h AND RETURN THE INCREMENTED VALUES AS yout, WHICH NEED NOT BE
  27. C A DISTINCT ARRAY FROM y.
  28. C SUBROUTINE derivs(x,y,dydx) RETURNS DERIVATIVES dydx AT x FOR n
  29. C DIFFERENTIAL EQUATIONS
  30. C -------------------------------------------------------------------------
  31.       
  32.       subroutine rk4 (y,dydx,n,x,h,yout)
  33. C     ==============  
  34.  
  35.       implicit real*4 (a-h,o-z)
  36.       parameter (nmax=10) 
  37.       dimension y(n),dydx(n),yout(n),yt(nmax),dyt(nmax),dym(nmax) 
  38.  
  39.       hh=h*0.5
  40.       h6=h/6. 
  41.       xh=x+hh 
  42.  
  43.       do 11 i=1,n 
  44.         yt(i)=y(i)+hh*dydx(i) 
  45. 11    continue
  46.  
  47.       call lorenz(xh,yt,dyt)
  48.  
  49.       do 12 i=1,n 
  50.         yt(i)=y(i)+hh*dyt(i)
  51. 12    continue
  52.  
  53.       call lorenz(xh,yt,dym)
  54.  
  55.       do 13 i=1,n 
  56.         yt(i)=y(i)+h*dym(i) 
  57.         dym(i)=dyt(i)+dym(i)
  58. 13    continue
  59.  
  60.       call lorenz(x+h,yt,dyt) 
  61.  
  62.       do 14 i=1,n 
  63.         yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
  64. 14    continue
  65.  
  66.       return
  67.       end 
  68.  
  69.  
  70. C -------------------------------------------------------------------------
  71. C RUNGE-KUTTA QUALITY CONTROL
  72. C FIFTH-ORDER RUNGE-KUTTA STEP WITH MONITORING OF LOCAL TRUNCATION ERROR
  73. C TO ENSURE ACCURACY AND ADJUST STEPSIZE. INPUT ARE THE DEPENDENT VARIABLE
  74. C VECTOR y OF LENGTH n AND ITS DERIVATIVE dydx AT THE STARTING VALUE OF THE
  75. C INDEPENDENT VARIABLE x. ALSO INPUT ARE THE STEPSIZE TO BE ATTEMPTED htry
  76. C THE REQUIRED ACCURACY eps, AND THE VECTOR yscal AGAINST WHICH THE ERROR
  77. C IS SCALED. ON OUTPUT, y AND x ARE REPLACED BY THEIR NEW VALUES, hdid IS 
  78. C THE STEPSIZE WHICH WAS ACTUALLY ACCOMPLISHED, AND hnext IS THE ESTIMATED
  79. C NEXT STEPSIZE. 
  80. C -------------------------------------------------------------------------
  81.  
  82.       subroutine rkqc (y,dydx,n,x,htry,eps,yscal,hdid,hnext)
  83. C     ===============
  84.  
  85.       implicit real*4 (a-h,o-z)
  86.       parameter (nmax=10,fcor=.0666666667,one=1.,safety=0.9,
  87.      &           errcon=6.e-4) 
  88.  
  89. C errcon EQUALS  (4/safety)**(1/pgrow)
  90.  
  91.       dimension y(n),dydx(n),yscal(n),ytemp(nmax),ysav(nmax),dysav(nmax)
  92.  
  93.       pgrow=-0.20 
  94.       pshrnk=-0.25
  95.       xsav=x
  96.  
  97.       do 11 i=1,n 
  98.         ysav(i)=y(i)
  99.         dysav(i)=dydx(i)
  100. 11    continue
  101.  
  102.       h=htry
  103. 1     hh=0.5*h
  104.  
  105.       call rk4(ysav,dysav,n,xsav,hh,ytemp) 
  106.  
  107.       x=xsav+hh 
  108.       call lorenz(x,ytemp,dydx) 
  109.       call rk4(ytemp,dydx,n,x,hh,y)
  110.  
  111.       x=xsav+h
  112. C      if(x.eq.xsav)pause 'stepsize not significant in rkqc.'
  113.  
  114.       call rk4(ysav,dysav,n,xsav,h,ytemp)
  115.  
  116.       errmax=0. 
  117.       do 12 i=1,n 
  118.         ytemp(i)=y(i)-ytemp(i)
  119.         errmax=max(errmax,abs(ytemp(i)/yscal(i))) 
  120. 12    continue
  121.       errmax=errmax/eps 
  122.  
  123.       if(errmax.gt.one) then
  124.         h=safety*h*(errmax**pshrnk) 
  125.         goto 1
  126.       else
  127.         hdid=h
  128.         if(errmax.gt.errcon)then
  129.           hnext=safety*h*(errmax**pgrow)
  130.         else
  131.           hnext=4.*h
  132.         endif 
  133.       endif 
  134.  
  135.       do 13 i=1,n 
  136.         y(i)=y(i)+ytemp(i)*fcor 
  137. 13    continue
  138.  
  139.       return
  140.       end 
  141.  
  142.  
  143. C -------------------------------------------------------------------------
  144. C RUNGE-KUTTA DRIVER WITH ADAPTIVE STEPSIZE CONTROL.
  145. C INTEGRATE THE nvar STARTING VALUES ystart FROM xstart WITH ACCURACY
  146. C eps TAKING nstep STEPS. 
  147. C hdes IS THE DESIRED STEPSIZE FOR WHICH THE RESULTS ARE WRITTEN TO FILE,
  148. C hmin IS THE MINIMUM ALLOWED STEPSIZE (CAN BE ZERO). 
  149. C -------------------------------------------------------------------------
  150.  
  151.       subroutine odeint 
  152.      & (ystart,nvar,xstart,nskip,nstep,path,eps,hdes,hmin, psig, pb, pr) 
  153. C     =================
  154.  
  155.       implicit real*4 (a-h,o-z)
  156.  
  157.       parameter (nmax=10,maxstp=10000)
  158.       parameter (two=2.0,zero=0.0,tiny=1.e-30) 
  159.  
  160.       dimension ystart(nvar),yscal(nmax),y(nmax),dydx(nmax)
  161.       dimension path(nmax,maxstp)
  162.  
  163.       common /parm/ sig,b,r
  164.  
  165.       x=xstart
  166.       kount=0
  167.  
  168.     sig = psig
  169.     b = pb
  170.     r = pr
  171.  
  172.       do 11 i=1,nvar
  173.         y(i)=ystart(i)
  174. 11    continue
  175.  
  176. C LOOP OVER STEPS
  177.       do 16 nstp=1,nskip+nstep
  178.  
  179.         h=hdes
  180.         hsum=0.
  181. 20      continue
  182.  
  183. C SET SCALE USED TO MONITOR ACCURACY
  184.         call lorenz(x,y,dydx) 
  185.         do 13 i=1,nvar
  186.           yscal(i)=abs(y(i))+abs(h*dydx(i))+tiny
  187. 13      continue
  188.  
  189. C CHECK COUNTER
  190.         kount=kount+1
  191. C        if(kount.gt.maxstp) pause 'max. number of steps exceeded.'
  192.  
  193. C GET NEW VALUES
  194.         call rkqc(y,dydx,nvar,x,h,eps,yscal,hdid,hnext)
  195. C        if(abs(hnext).lt.hmin) pause 'stepsize smaller than minimum.' 
  196.         hsum=hsum+hdid
  197.  
  198. C ONE MORE STEP ?
  199.         if(hsum.lt.hdes)then
  200.           h=min((hdes-hsum),hnext)
  201.           goto 20
  202.         endif      
  203.  
  204. C STORE INTERMEDIATE RESULT
  205.         if ( nstp.gt.nskip ) then
  206.           do 15 i=1,nvar
  207.             path(i,nstp-nskip)=y(i)
  208. 15        continue
  209.         endif
  210.   
  211. 16    continue
  212.  
  213.       return
  214.       end 
  215.  
  216.